{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "from math import factorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1) Polynomial Approximation with Derivatives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Brief Intro to `sympy`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we import `sympy`, a package for symbolic computation with Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we make a (symbolic) variable $x$ from which we can then build more complicated expressions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.var(\"x\")\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build up an expression with $x$. Assign it to `expr`. Observe that this expression isn't evaluated--the result of this computation is some Python data that represents the computation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sp.sin(sp.sqrt(x)+2)**2\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, take a derivative, using `.diff(x)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g.diff(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g.diff(x, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `.subs(x, ...)` and `.evalf()` to evaluate your expression for $x=1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g.subs(x,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g.subs(x, 1).evalf()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Helper function that takes symbolic functions as argument and plot them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_sympy(my_f, my_pts, **kwargs):\n",
    "    f_values = np.array([my_f.subs(x, pt) for pt in my_pts])\n",
    "    plt.plot(my_pts, f_values, **kwargs)\n",
    "\n",
    "    \n",
    "    \n",
    "def semilogy_sympy(my_f, my_pts, **kwargs):\n",
    "    f_values = np.array([my_f.subs(x, pt) for pt in my_pts])\n",
    "    plt.semilogy(my_pts, f_values, **kwargs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = 1/(20*x-10)\n",
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f.diff(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f.diff(x,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write out the degree 2 Taylor polynomial about 0:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "taylor2 = (\n",
    "    f.subs(x, 0)\n",
    "    + f.diff(x).subs(x, 0) * x\n",
    "    + f.diff(x, 2).subs(x, 0)/2 * x**2\n",
    ")\n",
    "taylor2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the exact function `f` and the taylor approximation `taylor2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pts = np.linspace(-0.4,0.4)\n",
    "\n",
    "plot_sympy(taylor2, pts, label=\"taylor2\")\n",
    "plot_sympy(f, pts, label=\"f\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.axis('equal')\n",
    "plt.grid()\n",
    "plt.xlabel('$x$')\n",
    "plt.ylabel('function values')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Behavior of the Error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Let's write the taylor approximation for any degree `n`, and define the error as f - tn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 2\n",
    "\n",
    "tn = 0\n",
    "for i in range(n+1):\n",
    "    tn += f.diff(x, i).subs(x, 0)/factorial(i) * x**i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_sympy(tn, pts, label=\"taylor\")\n",
    "plot_sympy(f, pts, label=\"f\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.axis('equal')\n",
    "plt.grid()\n",
    "plt.xlabel('$x$')\n",
    "plt.ylabel('function values')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We cand define the error as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "error = f - tn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_sympy(error, pts, label=\"error\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.ylim([-1.3, 1.3])\n",
    "plt.axis('equal')\n",
    "plt.grid()\n",
    "plt.xlabel('$x$')\n",
    "plt.ylabel('error')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get a better idea of what happens close to the center:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot only points close to zero [10^(-3),10^(0.1)]\n",
    "plt.figure(figsize=(10,6))\n",
    "pos_pts = 10**np.linspace(-3, 0.1)\n",
    "err_values = [error.subs(x, pt) for pt in pos_pts]\n",
    "plt.plot(pos_pts, err_values,'o')\n",
    "plot_sympy(tn, pos_pts, label=\"f\")\n",
    "plot_sympy(f, pos_pts, label=\"f\")\n",
    "plt.grid()\n",
    "plt.xlabel(\"$x$\")\n",
    "plt.ylabel(\"Error\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do you know why the discontinuity for the exact function? Take a look back at the analytical function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can take a look at the error now using loglog plot. We need to make sure we take the absolute value of the error first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.loglog(pos_pts, np.abs(err_values))\n",
    "plt.grid()\n",
    "plt.xlabel(\"$x$\")\n",
    "plt.ylabel(\"Error\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the slope of the error plot? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = sp.sqrt(x-10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 3\n",
    "x0 = 12\n",
    "\n",
    "tn = 0\n",
    "for i in range(n+1):\n",
    "    tn += f.diff(x, i).subs(x, x0)/factorial(i) * (x-x0)**i\n",
    "tn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The error of the Taylor approximation of degree 3 about x0 = 12 when h=0.5 is (that is, x = 12.5):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f.subs(x, 12.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tn.subs(x, 12.5).evalf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "error1 = f.subs(x, 12.5) - tn.subs(x, 12.5).evalf()\n",
    "abs(error1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now predict the error at $12.25$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and the actual error is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "error2 = f.subs(x, 12.25) - tn.subs(x, 12.25).evalf()\n",
    "abs(error2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 3:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sp.exp(x)\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Get the Taylor Series about the point $x_0 = 2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 4\n",
    "xo = 2\n",
    "\n",
    "taylor = 0\n",
    "\n",
    "for i in range(n):\n",
    "    taylor += g.diff(x, i).subs(x, xo)/factorial(i) * (x-xo)**i\n",
    "\n",
    "error =  g - taylor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pts = np.linspace(-1, 4, 100)\n",
    "plot_sympy(taylor, pts, label=\"taylor n=3\")\n",
    "plot_sympy(g, pts, label=\"f\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.grid()\n",
    "plt.xlabel('$x$')\n",
    "plt.ylabel('function values')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "semilogy_sympy(error, pts, label=\"error\")\n",
    "f2=x**2\n",
    "f3=x**3\n",
    "f4=x**4\n",
    "f5=x**5\n",
    "semilogy_sympy(f2, pts, label=\"$x^2$\")\n",
    "plot_sympy(f3, pts, label=\"$x^3$\")\n",
    "plot_sympy(f4, pts, label=\"$x^4$\")\n",
    "plot_sympy(f5, pts, label=\"$x^5$\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.grid()\n",
    "plt.xlabel('$x$')\n",
    "plt.ylabel('error')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "semilogy_sympy(error, pts, label=\"error\")\n",
    "f2=abs((x-2)**2)\n",
    "f3=abs((x-2)**3)\n",
    "f4=abs((x-2)**4)\n",
    "f5=abs((x-2)**5)\n",
    "semilogy_sympy(f2, pts, label=\"$(x-2)$\")\n",
    "semilogy_sympy(f3, pts, label=\"$(x-2)^3$\")\n",
    "semilogy_sympy(f4, pts, label=\"$(x-2)^4$\")\n",
    "semilogy_sympy(f5, pts, label=\"$(x-2)^5$\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.grid()\n",
    "plt.xlabel('$x$')\n",
    "plt.ylabel('error')\n",
    "plt.xlim([1.5,2.5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2) (Forward) Finite Difference Method "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return np.exp(x) - 2\n",
    "\n",
    "def df(x):\n",
    "    return np.exp(x)\n",
    "\n",
    "def df2(x):\n",
    "    return np.exp(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1, 2, 100)\n",
    "plt.plot(x, f(x), lw=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's evaluate the first derivative using finite difference approximation for decreasing values of h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xx = 1.0\n",
    "h = 1.0\n",
    "errors = []\n",
    "hs = []\n",
    "\n",
    "dfexact = df(xx) \n",
    "\n",
    "\n",
    "for i in range(20):\n",
    "    dfapprox = (f(xx+h) - f(xx)) / h\n",
    "    \n",
    "    err = np.abs(dfexact - dfapprox)\n",
    "    \n",
    "    print(\" %E \\t %E \" %(h, err) )\n",
    "    \n",
    "    hs.append(h)\n",
    "    errors.append(err)\n",
    "    \n",
    "    h = h / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.loglog(hs, errors, lw=3)\n",
    "plt.xlabel('h')\n",
    "plt.ylabel('error')\n",
    "plt.xlim([1e-6,1])\n",
    "plt.ylim([1e-6,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happens if we keep decreasing the perturbation h?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xx = 1.0\n",
    "h = 1.0\n",
    "errors = []\n",
    "hs = []\n",
    "\n",
    "dfexact = df(xx) \n",
    "fxx = f(xx)\n",
    "print('f exact = ',fxx)\n",
    "\n",
    "for i in range(60):\n",
    "    \n",
    "    fxxh = f(xx+h)\n",
    "    \n",
    "    dfapprox = (fxxh - fxx) / h\n",
    "    \n",
    "    err = np.abs(dfexact - dfapprox)   \n",
    "    \n",
    "    print(\" %E \\t %E\\t %E\" %(h, fxxh-fxx, err) )\n",
    "    hs.append(h)\n",
    "    errors.append(err)\n",
    "    \n",
    "    h = h / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.loglog(hs, errors, lw=3)\n",
    "plt.xlabel('h')\n",
    "plt.ylabel('error')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.loglog(hs, errors, lw=3,label='total')\n",
    "plt.loglog(hs,np.array(hs)*0.5*np.exp(1),'--',label='truncation')\n",
    "plt.loglog(hs,2*2.2e-16/np.array(hs),'--',label='rounding')\n",
    "plt.legend(bbox_to_anchor=(1.1, 1.05))\n",
    "plt.xlabel('h')\n",
    "plt.ylabel('error')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
